function Lshape
clc
% close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
projdir = '/home/hbui/kratos_janosch';

projname = 'Lshape.gid';
fn = [projdir '/' projname '/' 'mesh1.msh'];
% du = -0.01; lambda = [7.5 7.6:0.01:8];
du = -0.01; lambda = [7.5 7.6:0.05:25];

mg = fem2.gid.GidMeshGenerator(fn);
% mg = fem2.gmshmeshgenerator(fn);

mesh = mg.getMesh(2);
mesh.setNumberOfFields(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
meshutil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]);
c2 = fem2.Constraint([1 0]);
c2.setValue(1,du);
mesh.setConstraint(c1,meshutil.seekNodesOnSurface(0,1,-1000));
meshutil.seekNode([1000 500]).setConstraint(c2);
% plotMesh(mesh,{'trimesh','conditions'});
% return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = fem2.Model(mesh);
nf = mesh.countFaces();
%plain concrete data
e = 25850;
nu = 0.18;
t = 200;
ft = 2.7;
Gf = 0.065;
ef = Gf/ft;
Dnn = e/10;
Dss = e/10;
C = C_plane_stress(e,nu);
% mm = fem2.material.StVenantKirchhoffMM(0,0,t,0.0,fem2.enu.State.PLANE_STRESS);
for i=1:nf
    e = fem2.MatlabStructuralElement(mesh.getFace(i-1));
%     e = fem2.MatlabBBarStructuralElement(mesh.getFace(i-1));
%     e.setMaterialModel(mm);
    m.addElement(e);
    %create struct to contain element state variables
    es = createElementContainer();
    es.thickness = t;
    es.Dnn = Dnn;
    es.Dss = Dss;
    es.constitutiveMatrix = C;
    es.ef = ef;
    es.ft = ft;
    elements(i) = es;
end
%deactivate and activate elements
% for i=1:nf
%     c = computeFaceCenter(mesh.getface(i-1));
%     if (c(1)>=0 && c(1)<=0.6 && c(2)>=-0.6 && c(2)<=-0.4)
%         elements(i).isactivated = 1;
%     else
%         elements(i).isactivated = 0;
%     end
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
postprocessinfo.infoList = {'check sigma1','crack opening linder'};
postprocessinfo.infoParam = {[],[]};

% postprocessinfo.infolist = {'crack opening linder'};
% postprocessinfo.infoparam = {[]};

while 1
    [elements,output] = solving_strategy_static_2(elements,m,...
        lambda,50,1e-6,postprocessinfo);
    
    % [elements,output] = solving_strategy_static_4(elements,m,...
    %     lambda,dlambda,nt,50,1e-6,postprocessinfo);
    
    plotMeshWithCrack(m,elements,{'trimesh','conditions'});
    % plotdeformedmesh(mesh,0.2);
    
    disp('do you want to continue the analysis? ');
    conti = input('y(1)/n(0):');
    if conti == 1
        disp('analysis will continue from last lambda. please provide new lambda target:');
        lambda_max = input('next lambda target = ');
        last_lambda = lambda(size(lambda(:),1));
        step = input('lambda step = ');
        lambda = (last_lambda+step):step:lambda_max;
    else
        break;
    end
end
rmpath('../');
rmpath('../../');
